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Abstract. We report quantitative analysis of the radial gradient of solar angular 
velocity at depths down to about 15 Mm below the solar surface for latitudes up 
to 75° using the Michelson Doppler Imager (MDI) observations of surface gravity 
waves (f modes) from the Solar and Heliospheric Observatory (SoHO). A negative 
outward gradient of around — 400nHz/i?o, equivalent to logarithmic gradient of the 
rotation frequency with respect to radius which is very close to — 1, is found to be 
remarkably constant between the equator and 30° latitude. Above 30° it decreases 
in absolute magnitude to a very small value at around 50°. At higher latitudes the 
gradient may reverse its sign: if so this reversal takes place in a thin layer extending 
only 5 Mm beneath the visible surface, as evidenced by the most superficial modes 
(with degrees I > 250). The signature of the torsional oscillations is seen in this 
layer, but no other significant temporal variations of the gradient and value of the 
rotation rate there are found. 



1. Introduction 

The velocity field of the rotational flow in the Sun's near-surface layers 
may play a significant role in small-scale dynamo action in that region 
and in the dynamics of supergranular convection. Surface observations 
over decades and even centuries have shown that the latitudinal varia- 
tion of the surface rotation is rather smooth, being rather well described 
by a three-term (i.e. second-order) polynomial in fi 2 where \i = cos 9 
and 9 is the colatitude. Recent analyses of high-resolution data, in 
particular those utilizing solar f-mode observations by the Michelson 
Doppler Imager (MDI) on board the Solar and Heliospheric Obser- 
vatory (SoHO), have highlighted important departures from such a 
description of the rotation of the near-surface layers. The polar subsur- 
face layers (i.e. 6 < 20° and depths down to 28Mm below the surface) 
have been shown to be approximately lOnHz slower than expected 
from a simple three-term extrapolation from lower latitudes (Birch and 
Kosovichev, 1998; Schou et al., 1998; Schou, 1999); and Kosovichev and 
Schou (1997) have shown that, at a depth of 2 to 9 Mm beneath the 
surface, there exist zonal bands of alternate faster and slower rotation 
rate of ~ ±5m/s superimposed on the general trend described by the 
second order polynomial. This latter feature, inferred from the first 
observations of MDI in 1996, was found to be similar to the surface 
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'torsional oscillations' (Howard and Labonte, 1980) also observed in 
1995 in Doppler measurements using the first GONG observations 
(Hathaway et al., 1996).) More recently still, analysis of both p and f 
modes from the GONG network and MDI instrument have further led 
to the conclusion that these banded structures extend at least down to 
60Mm below the surface (Howe et al., 2000). 

The observed f modes, being confined to the outer layers of the Sun, 
provide a relatively clean and straightforward measure of conditions 
there. But those results above that were obtained just from the f modes 
assumed at least implicitly that the angular velocity is not varying 
significantly with depth within the layer sensed by those modes. It is 
however well known that another important property of the subsurface 
layers is that they present a radial gradient of angular velocity. This 
was first suggested by the fact that different indicators such as Doppler 
shifts of photospheric Fraunhofer lines, various magnetic field features 
of different ages and sizes (sunspots, faculae, network elements, H a 
filaments) or the supergranular network, present different rotation rates 
(see the review of Howard, 1984; Schroeter, 1985; Snodgrass, 1992). 
This has been interpreted by assuming that the different magnetic 
features are anchored at different depths (e.g. Foukal, 1972; Collin 
et al., 1995), their different rotation rate being therefore interpreted as 
an indication of the existence of radial gradients of angular velocity in 
the subsurface layers. More specifically, noticing that the supergranular 
network rotation rate (~ 473nHz) was found to be ~ 4% faster than the 
upper photospheric plasma rate obtained from spectroscopic methods 
and also ~ 2% faster than various magnetic indicators thought to be 
rooted under the supergranulation layer, Snodgrass and Ulrich (1990) 
inferred that a maximum of angular velocity should exist somewhere 
between O.95i? and the surface. 

From the theoretical point of view, it has been suggested that the 
angular momentum per unit mass Qr 2 sin 2 6 could be conserved in the 
supergranular flow (Foukal and Jokipii, 1975; Foukal, 1977; Gilman and 
Foukal, 1979). From dVl/VL = —2dr/r, at fixed latitude, this simple 
argument leads effectively to a negative gradient below the surface, 
and the 4% difference in rotation rates would be explained if the su- 
pergranulation network velocity observed at the surface were reflecting 
the rotation rate at a depth of 2%R Q ~ 15Mm, which turns out to 
correspond to the depth expected for the supergranular convection 
(Foukal, 1977; Duvall, 1980) (but see also Beck and Schou (2000) for 
more recent estimate). In order to reproduce the observed patterns of 
solar activity such as the equatorward migration of sunspots, early 
dynamo models based on a positive surface a-effect indicated also 
that the angular velocity must decrease outwards i.e. dQ/dr < (e.g. 
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Leighton, 1969; Roberts and Stix, 1972). One of the first goals of he- 
lioseismology was therefore to test the assumptions about the negative 
gradient of angular velocity below the surface suspected from differ- 
ent surface observations. This was indeed first attempted by Deubner 
et al. (1979): although they did not resolve individual modes, they 
were able, from ridge-fitting separately the eastward- and westward- 
propagating near-equatorial waves in the (k, w) diagram, to detect such 
a negative gradient close to the surface. Subsequent helioseismic work 
using resolved mode frequencies has shifted much theoretical focus to 
the base of the convection zone by showing that the radial gradient 
of angular velocity in the bulk of the convection zone is weak and 
that a strong radial shear, the so-called tachocline, occurs at its base. 
The gradient dtt/dr is positive in the tachocline at sunspot (i.e. low) 
latitudes (Brown et al., 1989). This has led various dynamo theories to 
locate the dynamo action below the convection zone, with a negative 
a-effect operating there (e.g. Gilman et al., 1989; Parker, 1993) though 
some recent work has revisited the idea of a positive surface a-effect 
but invoking the action of a meridional circulation, equatorward be- 
low the convection zone and poleward at the surface, to produce the 
observed equatorward migration of sunspots by advective transport of 
flux (Dikpati and Charbonneau, 1999; Kiiker et al., 2001). The lack 
until recently of precise determinations of high-degree mode param- 
eters made it difficult to obtain very localized inferences about the 
subsurface layers. But, because all the observed modes have amplitude 
close to the surface, inverters again got hints about the existence of 
a radial shear close to the surface (especially using methods such as 
regularized least-squares which readily extrapolate into regions where 
the data provide no localized information) though without being able 
to quantify precisely its extent and amplitude (e.g. Thompson et al., 
1996; Corbard et al., 1997). 

We show in this work that f-mode observations allow us to make 
quantitative inferences about the surface radial shear. These should be 
taken into account when modeling near-surface dynamo action or the 
dynamics of the supergranulation layer. 



2. Observations 

The data used here are 23 independent times series of 72 days obtained 
from the so-called MDI medium-/ program. These cover the period 
from 1996 May 1 to 2001 April 4 with interruptions during the sum- 
mer 1998 (June 23 to October 23) and between 1998 December 4 and 
1999 February 4 due to SoHO spacecraft problems. More details on the 
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production of these time series from the observations can be found in 
Schou (1999). 

A given f-mode multiplet in the spectra comprises 21 + 1 frequen- 
cies ui m , where I and m are the degree and azimuthal order of the 
spherical harmonic YJ m (#, 4>) describing the angular dependence of the 
modes. The so-called a coefficients for the multiplet are defined by the 
polynomial expansion: 

21 

vim = m + Y. al j V ?( m ) m = ±l,±2...±l (1) 
i=i 

where V are orthogonal polynomials normalized such that Vj l \l) = I 
(Schou et al., 1994). All f modes considered here have degrees between 
/ = 117 and / = 300 but the total number of multiplets observed is 
between 112 and 143, depending on the 72-day interval considered. 
For each observed mode, the central frequencies f/o and the first 36 a 
coefficients have been estimated using the method described in Schou 
(1992). Odd- indexed a coefficients, which describe the dependence of 
the frequencies that is an odd function of m, arise from the North- 
South symmetric part of the solar rotation. Even-indexed coefficients 
arise from latitudinal structural variation, centrifugal distortion and 
magnetic fields. 



3. Data analysis 

Following Ritzwoller and Lavely (1991), we identify the North-South 
symmetric part of the angular velocity f2(r, /i) with the odd-degree, 
zonal part of the toroidal component of a general stationary and lami- 
nar velocity field and write 

oo 

n(r^)=Y / ^ + i(r)T^), (2) 

j=0 

where r is fractional radius and T^- =T 2 1 J (/z)/T 2 1 ? (0) are Gegenbauer 
polynomials (see Appendix) normalized such that the equatorial rate 
is given by the straight sum of the f^j+iM- 

Assuming slow rotation, we can use a linear perturbation theory to 
predict the effect of rotation on the oscillation modes (e.g. Hansen et al., 
1977). Moreover, with the polynomials V and expansion Equation (2) 
as chosen, there is a one-to-one relation between odd a coefficients 
and the components ^j+iir) (Ritzwoller and Lavely, 1991), thereby 
reducing the full 2D problem to a set of ID integral equations often 
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refered as the 1.5D problem. In the particular case of the f modes, we 
obtain 



27r<4 +1 = u l 2j+1 J K l h (r)n 2j+ i(r)dr , 



(3) 



where the expression for the kernels K l h (r) and u l 2 j +1 are derived in 
Appendix. 
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Figure 1. f-modes rotational kernels K l h {r) for Z=117, 150, 200, 250, 300 from left 
to right. 

The 36 a coefficients extracted from observation do not provide infor- 
mation about the terms above j = 17 in the summation in Equation (2) 
and that corresponds to a limitation in the latitudinal resolution we can 
reach. Defining 



°2j+l 



2j+l 



W. 



I ' 
2j+l 

from Equations (2), (3) and (4) we obtain 

17 _ x 

E&W^M^/ K l h (r)n(r,fi )dr, 

3=0 ° 



(4) 



(5) 



where £l(r, /x) refers to the part of the rotation profile corresponding to 
the sum Equation (2) truncated at j = 17. We can also show (Pijpers, 
1997) that the above linear combination of b coefficients is such that 



17 _ /"I /"I 

J2 b 2j+i T 2j(^o) = / K l h (r)K(no,fj,)n(r,n)drdn , 

j—Q JO JO 



(6) 



where K(no,fi) are the so-called latitudinal averaging kernels which 
show what latitudinal average of the true rotation rate is made at 
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Figure 2. Latitudinal averaging kernels at 0, 20, 40, 60, 80° of latitude (dou- 
ble-dot-dash, dot-dash, dash, dot and full lines respectively) corresponding to the 
combination a Equation (8), b Equation (9). 



each latitude. Figure. 2a shows that these kernels have their main 
peak centered at ^lq but present an oscillatory behaviour which may 
lead to systematic errors if some small-scale features (corresponding to 
terms with j > 17) exist in the true rotation rate. In order to avoid 
this, one may try to find instead the combination of b coefficients that 
leads to kernels that are optimally localized around a given latitude. 
This can be achieved following for instance the method of Backus and 
Gilbert (1968), but we notice here that a similar result can be obtained 
simply by introducing, in the sum of Equation (6), a correcting factor 
e -i(i+3/2)// w h_ ere l Q = 117 corresponds to the lowest degree of the 
observed f modes (see also Equation (15)). Doing this, the latitudinal 
averaging kernels are found better peaked (Figure 2b) and the formal 
errors associated with the linear combination of the b coefficients is low- 
ered. Following the definition of Corbard et al. (2001), the latitudinal 
resolution obtained is about 10° at all latitudes. 

The kernels K l h (r) associated with each f mode have a simple shape 
with only one maximum located at slightly different radial positions 
depending on the degree / (Figure 1). If we define r l = Jq K l h (r)rdr, 
the radial location of the center of gravity of these kernels, and assume 
a linear behaviour of the rotation rate at each latitude in the radial 
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domain where the f modes considered have appreciable amplitude, i.e., 
ft(r,/i ) = a(/i )-/?(W))(r-l) (7) 
in r > 0.97, say, we simply obtain 

17 

E^+i^^o)^^,^), (8) 

3=0 

where the meaning of Cl is the same as in Equation (5). Alternatively, 
a slightly modified choice of weights yields 

17 

— JG/+3/2) 

Y, b 2j+i f 2j(^) e ^~ ~<^(^o^)>mo~ to(r l ,jio) , (9) 

3=0 

where the brackets denote the weighted average around the weight- 
ing function being the kernels of Figure 2b. The second approximate 
equality in Equation (9) would be exact if the rotation profile were a 
linear function of n 2 in the domain covered by the averaging kernels 
(i.e. ±10°), with /2q = Jq 1 k(/xo> /")/" 2 d/u; the approximation is less good, 
however, at high solar latitudes. 

The parameters a and (5 can then be estimated at each latitude from 
a linear least-square fit, yielding not only an estimate of the value of 
the rotation rate at e.g. the surface, but also an estimate of the average 
gradient dVL/dr in the region sampled by the f modes. Finally we note 
that the dependence of as a function of radius in the near-surface 
layers may sometimes conveniently be described by a power of r: we 
note that this description is immediately derivable from our a and /?, 
since for small values of 1 — r the right-hand side of Equation (7) is 
well approximated by a(fio)r~ a ^' 



4. Results 

By combining the frequency splittings within each f-mode multiplet 
in the manner given by Equation (8), for different choices of target 
latitude, we obtain measures of the near-surface rotation which are 
reasonably well localized in latitude and which correspond to different 
weightings in the depth direction. The latitudinal sensitivity is illus- 
trated in Figure 2 and the depth sensitivity in Figure 1. Figure 3 shows 
the results of combining the data using Equation (8), averaged in time 
over all the datasets under study. In depth, the points are plotted at the 
center of gravity (r = r l ) of the corresponding kernels (cf. Figure 1). It 
is evident from these results that, at low latitudes, the weighted rotation 
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Figure 3. Time average of f^rg, /io)/27r (Equation (8)) for values of /io corresponding 
to the latitudes indicated in each panel. The result of the linear fits (Equation (7)) are 
shown by the straight lines. The error bars are the standard deviation associated with 
the weighted temporal mean. The mark on the right of each panel indicate the surface 
plasma rate obtained by Snodgrass et al. (1984). Note that the surface spectroscopic 
value indicated on panel f is essentially an extrapolation from observations at lower 
latitudes. 

increases with depth. If at each latitude separately we fit these results 
to a rotation profile that is linear in depth, we obtain the linear fits 
overplotted in Figure 3. These provide an average rotational gradient 
(3(fio) in the outer 15 Mm or so of the solar interior, and an extrapolated 
surface rotation rate a(/io)- The gradient, as a function of latitude, is 
presented in Figure 4, both in terms of its dimensional value and in 
terms of the logarithmic derivative <91nf2/<91nr. It may be seen that 
for latitudes below 50° the gradient of rotation with depth is negative; 
at about 50° it is close to zero; and for higher latitudes the average 
rotational gradient becomes positive. We note that the radial gradient 
is remarkably constant at latitudes up to 30°, and the value of the 
logarithmic derivative at these latitudes is close to —1. We return to 
this in the Discussion. 

Another way to visualize the changing gradient with latitude is that 
in Figure 5, where we show the rotation rate extrapolated both to the 
surface (r = 1) and to r = 0.97. The deeper rotation is faster than the 
surface rotation at low- and mid-latitudes, but slower at high latitudes. 
At low- and mid-latitudes the extrapolated surface rate agrees well with 
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Figure 4- a Logarithmic derivative of angular velocity as a function of latitude. This 
corresponds to the ratio 13 /a of Equation (7). b Radial gradient of angular velocity 
f3 as a function of latitude, c Normalized \ value of the linear fit. The diamond 
symbols are for the results obtained using Equation (8) while the other points are 
obtained using Equation (9). The horizontal error bars indicate the angular resolu- 
tion as deduced from Figure 2b. The vertical error bars are formal errors deduced 
from the linear fit. The dashed horizontal line correspond to no radial gradient of 
angular velocity. 



the spectroscopic surface measurements, given the approximately 1.5% 
spread in recent such determinations (see the review by Beck (2000)). 
For comparison, we have made a fit to our inferred surface rate below 
60° latitude and present our fitting coefficients with the spectroscopic 
coefficients of Snodgrass et al. (1984) in Table I. Similarly to what has 
been found previously, our inferred rotation rate above 70° is markedly 
slower than what would be expected from a 3-term fit at low- and 
mid-latitudes: we return to this issue of the so-called 'slow pole' later. 

Figure 4c shows the chi-squared for the least-squares fits at each lat- 
itude. The large chi-squared values at higher latitudes are striking. The 
difference between the chi-squared values when using equations (8) and 
(9) is also very noticeable: this arises largely because the error bars on 
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Figure 5. The full line gives the photosheric plasma rotation rate inferred by Snod- 
grass et al. (1984); the diamond symbols and horizontal bars correspond to a, the 
intercept of the linear fit respectively in the case of Equation (8) and Equation (9) 
and the dashed line correspond to an extrapolation of the rotation rate at O.97i?0 
using Equation (7) in the case of Equation (8). 

the fitted points are reduced by the exponential factor in Equation (9), 
which results in an increased chi-squared. Thus the interpretation of 
the absolute value of the chi-squared may be a little uncertain, but the 
trend with latitude for the two cases is similar. The higher values of 
chi-squared at higher latitudes is consistent with the greater deviation 
from a linear fit in the high-latitude panels of Figure 3. The systematic 
deviation of the near-surface points contributes most to the chi-squared: 
these correspond to the high-degree modes and so motivates taking a 
closer look at those data. (The scatter of the deepest points is large 
but less significant because of the large error bars on those points.) 

We have therefore repeated our analysis but excluding those modes 
of degree / > 250 and / < 160. The resulting gradient and chi-squared 
are shown in Figure 6. Compared with the previous result (Figure 4) the 
gradient is similar for latitudes lower than 50 degrees. Now it is evident 
from Figure 3 that, at high latitudes, excluding the high-degree modes 
will tend to make the fitted gradient less positive. Indeed, we find that 
the gradient without the / > 250 data remains slightly negative up to 
about 75 degrees. Also, the values of chi-squared have been more than 
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Table I. Comparison between the surface plasma rate and our 
results from f modes analysis. 





fi!/27T 


ft 3 /27T 


n 5 /27r 


Method 


(nHz) 


(nHz) 


(nHz) 


Snodgrass et al. (1984) 1 


436.4 


21.0 


-3.6 


f modes (Z-averaged) 2 








(117 < I < 300) < r l >= 0.991 3 


438.8 


21.0 


-3.9 


(160 < I < 250) < r l >= 0.991 


438.9 


21.2 


-4.0 


f modes (surface extrapolation) 4 








(117 < I < 300) 


435.8 


20.2 


-3.2 


(160 < I < 250) 


435.7 


20.5 


-3.6 



1 Spectroscopic measurements made at the Mount Wilson 150- 
foot Tower between 1967 and 1982. 

2 Average of the first 3 b coefficients (cf. Equations (4), (8)). 
3 Center of gravity of the corresponding I averaged radial kernels. 
4 Obtained by fitting the intercept a(fi) to the expansion Equa- 
tion (2) for latitudes below 60°. 



halved at high latitude, compared with our previous linear fit to all the 
f-mode data (Figure 6b). The inferred low- and mid-latitude surface 
rate is barely affected (compare the last two lines of Table I) . 

It is interesting also to compare the linear fit to the I < 250 data 
with a fit of a constant function to the same data: the constant fit is 
equivalent averaging the f-mode splittings over I (see Table I). It is 
evident from Figure 6b (dotted line) that this provides a very poor 
fit below about 55°: the data strongly favour the model with a linear 
depth-dependence there. At high latitudes, the linear fit selects only 
a very small gradient and so the two chi-squared functions are very 
similar: the data for I < 250 indicate that at high latitudes the gradient 
is small, in the range of depths spanned by their lower turning points. 

If the data for I > 250 are indeed reliable, then the discrepancy 
between the results in Figures 4a and 6a implies that the model of 
rotation varying linearly with depth is not appropriate at high latitudes 
and the extrapolation to the surface at those latitudes will be unreli- 
able. An alternative approach is to attempt to construct kernels that 
are localized in depth using the Optimally Localized Averaging (OLA) 
kernel in depth (cf. Christensen-Dalsgaard et al., 1990) in the manner 
of (Backus and Gilbert, 1968). Such kernels at two selected depths 
are shown in Figure 7: they were constructed using all the available 
f modes. It should be noted that the method succeeds in producing 
kernels which are reasonably localized and which have their center of 
gravity outside the range of abscissa values in Figure 3, that is, the 
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Figure 6. Similar to Figure 4 but using only modes 160 < I < 250. The radial 
gradient of angular velocity remains negative even at latitudes above 55° . The dotted 
line on panel b shows, in the case of Equation (8), the \ 2 values corresponding to a 
fit by a constant which is equivalent of taking an average over I. 



method uses the mode sensitivities to extrapolate to greater depths and 
closer to the surface. In particular, in the latter case one expects that 
the increasing trend of values for the near-surface points in Figure 3e,f 
means that the near-surface Backus-Gilbert inversion at those latitudes 
will have values higher than those seen in Figure 3. This is exactly what 
is found (Figure 7): the Backus-Gilbert inversion at high latitudes for 
r = 0.986 interestingly falls below the 3-term spectroscopic surface 
rate, but even more strikingly the corresponding near-surface result at 
r = 0.997 lies above it by 2 — 4 standard deviations. This is another way 
of demonstrating that the increasing values of the combined splittings 
for I > 250, if they are reliable, indicate a strongly positive gradient of 
rotation with radius in the rather superficial subsurface layers at high 
latitudes. 

To look for possible temporal variations of the subsurface shear, we 
have analyzed each one of the 23 72-day datasets individually in exactly 
the same manner as we analyzed the time-averaged set (e.g. Figure 3), 
and derived an intercept value a(no;t) (corresponding to the surface 
rate at that location and epoch) and slope (3(no;t) from a linear fit 
to the combined splittings for each latitudinal location (1q and time 
t. The resulting estimated surface rates and slopes at three latitudes 
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Figure 7. Rotation profiles as a function of latitude corresponding to depth averaged 
shown in the sub panel. The dashed and full lines correspond respectively to the 
shallower and deeper kernels which have respectively O.986i? and O.997i? as center 
of gravity. The dot-dashed line corresponds to the Snodgrass et al. (1984) plasma 
rotation rate. These results are obtained by using all modes from I = 117 to I = 300. 

(equator, 30°, 60°) are shown in Figure 8. The large-scale variations 
in the surface rate correspond very well to the migrating banded zonal 
flows (torsional oscillations) measured by Schou (1999) and by Howe 
et al. (2000): the equatorial surface rate starts high because of the tail- 
end of one migrating band of faster flow, then drops down and rises 
again towards solar maximum as another band of faster flow reaches 
the equator: the latter was at 30° at the beginning of the cycle, hence 
the rate at that location starts high and drops as the band migrates 
closer to the equator. The 60° rate rises as the high-latitude banded 
flow reported by e.g. Schou (1999) strengthens towards solar maximum. 
The slope shows no significant corresponding variations, implying that 
the torsional oscillations raise and lower the rotation rate across the 
whole depth of the layer without changing the shear gradient. There 
are indications of annual variations in the inferred values of the slope 
(most strikingly at 30°), which are almost certainly an artifact: such 
artifacts can conceivably arise from annual variations in SoHO's orbit. 
Other evidence for one-year artifacts in the f-mode data is presented by 
Antia et al. (2001). These should not affect the time-averaged values, 
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Figure 8. Intercept (left column) and slope (right column) of the linear fit 
Equation (7) at the equator, 30° and 60° of latitude (from top to bottom). 



however. There is no noticeable annual variation in our inferred values 
of the surface rotation rate. 



5. Discussion 



We have used the depth and latitude variation in the sensitivities of 
the solar f modes to deduce the rotation profile in the subsurface shear 
layer of the Sun in the outer 15 Mm of the solar interior. Our work 
differs from earlier seismic investigations. These were either based on 
the f modes but implicitly assumed a depth-independent model of the 
rotation (e.g. Schou, 1999), or used global inversions of p- and f-mode 
splittings and consequently may suffer from any systematic difference 
between the p- and f-mode data (e.g. Schou et al., 1998), or used local 
helioseismic ring analysis (e.g. Basu et al., 1999; Haber et al., 2000), 
which promises to be a powerful diagnostic of near-surface flows and 
stratification but the sensitivity and systematics of which are still under 
investigation (Hindman et al. 2001, in preparation). By using just the 
splittings of the f modes, which are arguably the most straightforward 
helioseismic modes to interpret, we believe we are able to obtain not 
only a simple but also a clean measure of the near-surface shear. As 
with all inferences about rotation from global splittings, we note that 
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only that component of rotation which is symmetric about the equator 
is recovered. 

The most robust results concern the low-latitude shear. The aver- 
age gradient dlnQ/dlnr (at constant latitude) in the outer 15 Mm is 
close to —1 and remarkably constant from the equator to 30° latitude. 
Between 30° and 55° latitude, the gradient is still negative but makes 
a steady transition to a small (absolute) value. All our analyses show 
this. The variation of rotation at these latitudes appears to be well 
described by a linear function of depth, within the outer 15 Mm. 

As discussed in the Introduction, if moving parcels of fluid were to 
conserve their specific angular momentum as they moved towards or 
away from the rotation axis, one would find that the rotation rate varied 
as the inverse square of the distance from the axis of rotation, so at low 
latitudes one would have that d In 0,/d In r ~ —2. In reality other effects 
such as diffusion will cause exchange of angular momentum between 
parcels, so we may expect a logarithmic gradient somewhat smaller 
in magnitude than —2. A precise measurement of this value in the 
Sun provides information about the relative effectiveness of competing 
mechanisms transporting angular momentum. Our finding is that at 
latitudes below 30° the value of the logarithmic gradient is much closer 
to —1 than to —2. In fact, this seems in reasonable agreement with the 
equatorial value found by DeRosa in numerical simulations of rotat- 
ing compressible convective fluid in a thin shell representing the Sun 
between about 0.94R and 0.98.R (DeRosa, 2001; DeRosa et al., 2001). 
Also these simulations show a tendency for the gradient to decrease in 
magnitude as one moves from equator to mid-latitudes, albeit at lower 
latitudes than we find for the Sun. Although these simulations exclude 
for numerical reasons the near-surface layers that we are probing, the 
qualitative agreement is nonetheless encouraging. 

At latitudes above ~ 55°, the depth-averaged gradient over the layer 
appears to change sign with respect to the low-latitude shear, though 
this is largely a consequence of the behaviour in the very near-surface 
layers (outer 5 Mm) which in turn is deduced from the splittings of 
the highest-degree f modes. The gradient in the range of depths 5 — 
15 Mm is small at these high latitudes; and such significant gradient 
dfl/dr as does exist at high latitude (if any) is in the outer 5 Mm and 
predominantly positive. We note that, using a ring-analysis technique, 
Basu et al. (1999) deduced a similar behaviour at high latitudes, finding 
a reversal of gradient in a zone above O.994i?0. 

Concerning the surface rotation rate itself, below 55° our extrapola- 
tion of the rotation rate to the surface is in satisfactory agreement with 
the directly measured spectroscopic surface rotation rate (cf. Table I). 
Our inferred surface rate should be more accurate than one simply 
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inferred from the averaged f-mode splittings, because we take out the 
linear gradient with depth which undoubtedly exists at these latitudes: 
this can make a difference of ~ 5nHz, even over the fairly small range 
of depths sampled by the observed f modes. 

The seismically inferred surface rate at high latitudes is considerably 
less secure. It has previously been reported from helioseismic investi- 
gations that the high-latitude surface rate is lower what one would 
expect from a simple three-term extrapolation from lower latitudes 
(Schou et al., 1998; Schou, 1999). Indeed it can be seen from panel 
f of Figure 3 that many of the points fall below the extrapolated 
spectroscopic rate for that latitude, implying that the rotation rate 
at some depth is lower than the spectroscopic surface rate one would 
infer from the values in Table I. The rather flat plateau of values in 
those panels strongly suggests that the rotation rate at about 10-15 
Mm depth is slower than the extrapolated spectroscopic rate, which is 
confirmed by our OLA inversion result at those depths. However, the 
combined splittings at high degree are increasing with / and if taken at 
face value, as is done in our OLA inversion result for r = 0.997i?, this 
behaviour implies that the very near-surface rotation rate is actually 
higher than the spectroscopic rate. Thus the matter is still open. Since 
the quoted spectroscopic rate is principally an extrapolation of surface 
observations at low- and mid-latitudes, the true rotation rate that 
would be determined by spectroscopy at high latitudes is uncertain. 
Direct spectroscopic determinations at high latitude would resolve the 
question. The very high-degree splittings could contain some systematic 
errors, and if these affect the low-m data the most (some evidence for 
such an effect for p modes at lower degrees is offered by the comparison 
of GONG and MDI splittings by Schou et al. (2001)), then the near- 
surface, high-latitude rotation rates inferred here could be erroneously 
high. We hope that this possibility will shortly be addressed by inde- 
pendent determinations of these splittings by the GONG experiment 
using the new higher-resolution GONG+ observations. 



6. Conclusion 

Finally, to return again to our principal focus which is the shear gradi- 
ent of the near-surface rotation, we find that at low and mid-latitudes 
the gradient dVi/dr in the outer 15 Mm or so is close to —1 and is quite 
independent of latitude below 30°; between 30° and ~ 50° latitude, it is 
still negative but makes a transition to small absolute value. At higher 
latitudes, the gradient in the bulk of the outer 15 Mm is probably 
small, but if the highest-degree (/ > 250) data are to be believed there 
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is a region of positive gradient in the outer 5 Mm at high latitudes, 
similar to what Basu et al. (1999) found from ring analysis. We find 
no evidence for the gradient to vary with time: the torsional oscillation 
seems to pass through without changing the shear gradient in the outer 
15 Mm. 

Interestingly, the most recent circulation-dominated dynamo mod- 
els (Dikpati and Charbonneau, 1999; Kiiker et al., 2001) are able to 
reproduce to some extent the equatorward migration patterns without 
invoking any radial gradient of angular velocity at the surface. Such 
negative gradient at low latitude should however probably be taken 
into account because if it is associated with a positive surface a-effect, 
it will compete against the surface poleward circulation and contribute 
to producing the equatorward migration of magnetic patterns observed 
at the surface of the Sun. 
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Appendix 

Derivation of f-mode 1.5D kernels 

The polynomial V used to describe the frequency splittings can be 
expressed in terms of the Clebsh-Gordan coefficients C 3 ™ mi ^ 2rrl2 (e.g. 
Edmonds, 1960) by: 



The Gegenbauer polynomials used in Equation (2) are defined by (e.g. 
Morse and Feshbach, 1953): 



tUu)- ^ fin 
W-V47 + 3 • 
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From Ritzwoller and Lavely (1991) we can deduce that 



vL., r i 



27r °^ +1 = Tpj) Jo ^( r ) n 2i+i«dr , (12) 
where = L 2 C\\ {2j+l) J f3 l 2j+l , L 2 = 1(1 + 1) and 



£j and 77/ being respectively the radial and horizontal displacement 
eigenfunctions which are determined by solving the differential equa- 
tions describing the motion of a self-graviting fluid body in a standard 
solar model (e.g. Unno et al., 1989) and p is the density profile given 
by the model, all these being functions of the fractional solar radius r. 

Other expressions of practical interest can be found for w^'+i that are 
recalled here for completeness. Pijpers (1997) established the recurrence 
relation 

(j-Z)(2j + l) j 
j(2l + 2j + 1) 

and Schou (1999) noticed that to a very good approximation 



l u ~ VlfJ t ±i i . , 



4 J+1 /T 2 1 J .(0)«e-^^W. (15) 

The f modes are horizontally propagating surface gravity waves 
for which the displacement eigenfunctions satisfy the following surface 
boundary condition under the Cowling approximation (e.g. Berthomieu 
and Christensen-Dalsgaard, 1991): 

m(r) « ^6(r), (16) 

where g s = GM & /Rq is the surface gravitational acceleration. More- 
over, the angular frequencies wi = 2ttuiq of the f modes follow asymp- 
totically (for I — > oo) the dispersion relation wf « g s L/R & . Therefore 
we have & « and, from Equation (13), the rotational kernels asso- 
ciated with the f modes can be written as a function of the horizontal 
displacement only: 

K l h (r) = wSt)1pM r ' 2 



K l Ar) « k\K l h (r) { h 1 Jo^iWW rMr . (17) 

' ' fcj = l-i-^(l+i(2i + 3)) 



Finally, Equation (3) is obtained by taking: 



u l 2j+1 « fc{ e -^'+ 3 / 2 )/*. (18) 
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We note that Equation (5) is obtained by using the fact that, in the 
approximation Equation (17) valid for f modes, the rotational kernels 
depend on j only by a multiplicative factor. Taking instead K l - ~ Kq 
for all j as usually done for high degree modes would also allow us to 
write Equation (5) but the integrated difference j(Kj — K l )dr would 
reach 2.2% for / = 117, j = 17 whereas it remains negligible for all I 
and j in the case of the approximation used here. 
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